function [schidx,nsidx] = HiLo_TrialSelection(sch,schidx)
% [schidx,nsidx] = HiLo_TrialSelection(sch,schidx)
%
%  Use Gellermann's rules for two alternative choice task trial sequences
%  (Gllermann, 1933).

%  Takes the schedule and schidx as input.
%       SCH is the schedule structure
%       SCHIDX is a 1D matrix with the number of rows equal to the number
%           of unique trials; i.e. size(trials,1).
%  Returns an updated version of SCHIDX and NSIDX
%       NSIDX is the next schedule trial index, i.e. the next index
%           selected from the SCH.trials matrix to run.
%       SCHIDX should be updated with the number of runs for each unique
%           trial in SCH.trials.  The returned version should have
%           SCHIDX(NSIDX) incremented.
%
%  This same format for input and output paramters must be met for custom
%  trial select functions (i.e., experiment.trialselectfcn);
%
%  DJS (c) 2012

persistent GSEQ GSEQseed

global G_RP G_DA

if isempty(GSEQ)
    % Generate Gellermann (1933) randomized sequence of trials
    % Gellermann, J1933
    n = 10;
    p = zeros(2^14,n);
    for i = 1:n
        p(:,i) = randperm(2^14);
    end
    
    m = [-ones(2^12,n); ones(2^12,n)];
    m = unique(m(p),'rows'); clear p
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    % rule 1
    ind = sum(m,2) ~= 0;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    % rule 2
    tm1 = m > 0;
    tm2 = m < 0;
    for i = 1:n-4
        ind = sum(tm1(:,i:i+3),2) > 3;
        ind = ind | sum(tm2(:,i:i+3),2) > 3;
        tm1(ind,:) = [];
        tm2(ind,:) = [];
        m(ind,:)   = [];
    end
    
    %
    ind = false(size(m,1),n-2);
    tm1 = m > 0;
    tm2 = m < 0;
    for i = 1:n-2
        ind(:,i) = sum(tm1(:,i:i+2),2) == 3;
        ind(:,i) = ind(:,i) | sum(tm2(:,i:i+2),2) == 3;
    end
    ind = sum(ind,2) >= 2;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    % rule 3
    ind = abs(sum(m(:,1:n/2),2))>1;
    ind = ind | abs(sum(m(:,n/2+1:end),2))>1;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    % rule 4
    dm = diff(m,1,2)/2;
    ind = sum(abs(dm),2) ~= n/2;
    m(ind,:) = [];
%     imagesc(m); title(sprintf('n = %d',size(m,1)));
    
    % rule 5
    % rule 5 is implied
    
    % Construct trial sequence
    a = m(1:size(m,1)/2,:);
    b = m(size(m,1)/2+1:end,:);
    
    ridx = randperm(size(a,1)); a = a(ridx,:);
    ridx = randperm(size(b,1)); b = b(ridx,:);
    
    m = [a fliplr(b)];
    
    GSEQ = reshape(m',1,numel(m));

    % Start from random position within sequence
    % Since Gellermann rules for randomization only allow for 44 sequences
    % of 10 trials yield 440 trials total, limit seed to first 100 trials
    % in the sequence.
    GSEQseed = randperm(100);
    GSEQseed = GSEQseed(1);
end

if ~isfield(sch,'trial_idx'), sch.trial_idx = 1; end
trials = sch.trials;
boxid  = sch.boxid;

RV = sch.response_vals;

% correct side: -1 Left; 1 Right
csidx = findincell(strfind(sch.writeparams,'CorrSide')); 
corrside = cell2mat(sch.trials(:,csidx));


if isempty(RV)
    tidx = 0;
else
    tidx = size(RV,1);
end


subcorrside = find(corrside == GSEQ(GSEQseed + tidx));


% initialize first trial
if isempty(schidx)
    schidx = zeros(size(trials,1),1); % initialize SCHIDX
end


% give priority to least chosen trials
i = min(schidx(subcorrside));
i = find(schidx(subcorrside) == i);
r = randperm(length(i));
nsidx = subcorrside(i(r(1)));

schidx(nsidx) = schidx(nsidx) + 1;



% Specific for HiLo paradigm
pfidx = findincell(strfind(sch.writeparams,'ProbeFreq'));
rfidx = findincell(strfind(sch.writeparams,'RefFreq'));
ttidx = findincell(strfind(sch.writeparams,'trial_type'));

pfrange = trials{1,pfidx};
rfrange = trials{1,rfidx};

if trials{nsidx,ttidx} == 0 % No reference (Quiet trial)
    rf = 0;
    pf = randi(pfrange);
else
    rf = randi(rfrange);
    pf = rf * 2 ^ (trials{nsidx,csidx}*0.5); % half octave above/below reference
end

if isempty(G_DA)
    G_RP(sch.writemodule(rfidx)).SetTagVal(sprintf('RefFreq~%d',boxid),rf);
    G_RP(sch.writemodule(rfidx)).SetTagVal(sprintf('ProbeFreq~%d',boxid),pf);
else
    G_DA(sch.writemodule(rfidx)).SetTargetVal(sprintf('Stim.RefFreq~%d',boxid),rf);
    G_DA(sch.writemodule(rfidx)).SetTargetVal(sprintf('Stim.ProbeFreq~%d',boxid),pf);
end
    
% look for RewardRate parameter
rridx = findincell(strfind(sch.writeparams,'RewardRate'));
if ~isempty(rridx)
    if sum(schidx) > 10 % first 10 trials are always rewarded.
        r = rand <= sch.trials{nsidx,rridx}/100;
    else
        r = 1;
    end
    ststr = sprintf('RewardTrial~%d',boxid);
    if isempty(G_DA)
        G_RP(sch.writemodule(rridx)).SetTagVal(ststr,r);
    else
        G_DA(sch.writemodule(rridx)).SetTargetVal(['Stim.' ststr],r);
    end
end
